Evolution of clinical characteristics

1 Plot delta from baseline and evolution of significantly changed clincal parameters

xlabels <- c("0", "3", "6", "9", "12", "15", "18", "21", "24")

full_metadata <- as_data_frame(sample_data(ps_full))

full_metadata_distinct <- full_metadata%>% # to not have multiple data entries per visit/patient 
  distinct(id_visit, .keep_all = T)%>%
  filter(!is.na(visit_cal_cor))%>%
  mutate(baseline_ppfev1 = case_when(visit_cal_9=="1"~ pp_fev_percent))%>%
  group_by(id)%>%
  fill(baseline_ppfev1)%>%
  mutate(delta_ppfev1 = pp_fev_percent-baseline_ppfev1)%>%
  ungroup()%>%
  mutate(baseline_abtt_days = case_when(visit_cal_9 == "1" ~ treatmentdays_365prior_visit))%>%
  group_by(id)%>%
  fill(baseline_abtt_days)%>%
  mutate(delta_atb_tt_days = treatmentdays_365prior_visit-baseline_abtt_days)%>%
  ungroup() %>% 
  mutate(baseline_chloride = case_when(visit_cal_9 == "1" ~ chlorid_mmol_cl_l))%>%
  group_by(id)%>%
  fill(baseline_chloride)%>%
  mutate(delta_chloride = chlorid_mmol_cl_l-baseline_chloride)%>%
  ungroup() %>% 
  mutate(baseline_igG = case_when(visit_cal_9 == "1" ~ ig_g_g_l))%>%
  group_by(id)%>%
  fill(baseline_igG)%>%
  mutate(delta_igG = ig_g_g_l-baseline_igG)%>%
  ungroup()

#### delta ppfev1 ####
lm <- summary(lmerTest::lmer(delta_ppfev1~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_ppfev1)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Dppfev <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Dppfev $N_obs <- replace_na(lmm_Dppfev $N_obs, lmm_Dppfev $Total_N_obs[1])

lmm_Dppfev  <- lmm_Dppfev  %>% 
  mutate(Parameter = "Delta ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Dppfev1<- full_metadata_distinct%>%
    ggplot(aes(y=delta_ppfev1, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta ppFEV1")+
    xlab("Months from ETI treatment start")+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 43, step.increase = 0.05)
Dppfev1

#### delta chloride ####
lm <- summary(lmerTest::lmer(delta_chloride~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_chloride)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Dchloride <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Dchloride $N_obs <- replace_na(lmm_Dchloride $N_obs, lmm_Dchloride $Total_N_obs[1])

lmm_Dchloride  <- lmm_Dchloride %>% 
  mutate(Parameter = "Delta sweat chloride (mmol/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

lm_stats <- lm_stats %>% 
  filter(Months_after_ETI_start%in%c("3","12","18","24"))

Dchloride<- full_metadata_distinct%>%
  filter(visit_cal_9%in%c("1","2","5","7","9")) %>% 
    ggplot(aes(y=delta_chloride, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
        geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta sweat chloride (mmol/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= c("0","3","12","18","24"))+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 4, step.increase = 0.05)
Dchloride

#### delta IgG ####
lm <- summary(lmerTest::lmer(delta_igG~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_igG)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Digg <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Digg $N_obs <- replace_na(lmm_Digg $N_obs, lmm_Digg$Total_N_obs[1])

lmm_Digg  <- lmm_Digg %>% 
  mutate(Parameter = "Delta IgG (g/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Digg<- full_metadata_distinct%>%
    ggplot(aes(y=delta_igG, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
        geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta IgG (g/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.4, step.increase = 0.05)
Digg

#### delta antibiotic treatment days ####
lm <- summary(lmerTest::lmer(delta_atb_tt_days~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_atb_tt_days)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Datb <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Datb $N_obs <- replace_na(lmm_Datb $N_obs, lmm_Datb $Total_N_obs[1])

lmm_Datb  <- lmm_Datb%>% 
  mutate(Parameter = "Delta antibiotic tt days") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Datb<- full_metadata_distinct%>%
    ggplot(aes(y=delta_atb_tt_days, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta antibiotic tt days")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 23, step.increase = 0.05)
Datb

ggarrange(Dchloride, Dppfev1, Digg, Datb, nrow = 1)

lm_stats_1 <- rbind(lmm_Dppfev, lmm_Dchloride, lmm_Digg, lmm_Datb)

kable(lm_stats_1)

Parameter	Months_after_ETI_start	N_obs	Total_N_obs	Estimate	Std..Error	df	t.value	p	fdr	fdr_star
Delta ppFEV1	Baseline (Intercept)	35	235	1.0158422	3.4045717	38.64296	0.2983759	0.76701	0.84372	ns
Delta ppFEV1	3	30	235	12.4628627	1.5482261	193.34212	8.0497692	0.00000	0.00000	***
Delta ppFEV1	6	31	235	11.3119118	1.5384045	194.31464	7.3530151	0.00000	0.00000	***
Delta ppFEV1	9	31	235	13.4351520	1.5389186	194.47907	8.7302550	0.00000	0.00000	***
Delta ppFEV1	12	26	235	9.3934309	1.6321052	196.03302	5.7554077	0.00000	0.00000	***
Delta ppFEV1	15	22	235	10.2279562	1.7298648	196.79571	5.9125754	0.00000	0.00000	***
Delta ppFEV1	18	20	235	7.4812868	1.7930864	198.43422	4.1722956	0.00005	0.00008	***
Delta ppFEV1	21	22	235	6.8596951	1.7354919	198.80899	3.9525942	0.00011	0.00017	***
Delta ppFEV1	24	18	235	6.3404820	1.8748622	200.42886	3.3818389	0.00087	0.00119	**
Delta ppFEV1	sex	235	NA	-1.4902106	2.9866734	32.30319	-0.4989533	0.62119	0.75924	ns
Delta ppFEV1	age_y	235	NA	-0.0100933	0.1079618	33.07653	-0.0934894	0.92608	0.92608	ns
Delta sweat chloride (mmol/l)	Baseline (Intercept)	34	116	-0.4461111	4.0390632	40.67402	-0.1104492	0.91260	0.91260	ns
Delta sweat chloride (mmol/l)	3	24	116	-36.3784810	2.7750387	78.98507	-13.1091795	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	6	6	116	-34.7004282	4.9665859	90.79944	-6.9867771	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	9	3	116	-31.8304483	6.8780454	90.56285	-4.6278334	0.00001	0.00002	***
Delta sweat chloride (mmol/l)	12	18	116	-40.9264249	3.0860935	81.88236	-13.2615635	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	15	5	116	-38.3445596	5.3359683	88.97815	-7.1860546	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	18	14	116	-44.8826099	3.4289466	83.23005	-13.0893287	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	21	2	116	-31.6546927	8.1648932	87.76237	-3.8769267	0.00020	0.00028	***
Delta sweat chloride (mmol/l)	24	10	116	-42.4302600	3.9106492	82.87276	-10.8499274	0.00000	0.00000	***
Delta sweat chloride (mmol/l)	sex	116	NA	3.2813715	3.4512967	30.08667	0.9507648	0.34930	0.42692	ns
Delta sweat chloride (mmol/l)	age_y	116	NA	-0.0477268	0.1255767	31.08161	-0.3800613	0.70648	0.77713	ns
Delta IgG (g/l)	Baseline (Intercept)	30	200	1.0007796	0.7210072	32.36255	1.3880300	0.17461	0.19207	ns
Delta IgG (g/l)	3	26	200	-1.2117237	0.3043869	163.41537	-3.9808664	0.00010	0.00016	***
Delta IgG (g/l)	6	26	200	-1.5615271	0.3056767	164.34128	-5.1084264	0.00000	0.00000	***
Delta IgG (g/l)	9	26	200	-1.6885322	0.3060615	164.68197	-5.5169710	0.00000	0.00000	***
Delta IgG (g/l)	12	21	200	-1.7800419	0.3288477	165.79080	-5.4129674	0.00000	0.00000	***
Delta IgG (g/l)	15	18	200	-1.7938293	0.3478857	166.52434	-5.1563747	0.00000	0.00000	***
Delta IgG (g/l)	18	18	200	-1.5590856	0.3487345	167.71502	-4.4706945	0.00001	0.00003	***
Delta IgG (g/l)	21	18	200	-1.6565225	0.3486516	168.24383	-4.7512256	0.00000	0.00001	***
Delta IgG (g/l)	24	17	200	-1.2153766	0.3592677	170.20725	-3.3829275	0.00089	0.00122	**
Delta IgG (g/l)	sex	200	NA	0.3700112	0.6559188	27.74904	0.5641113	0.57721	0.57721	ns
Delta IgG (g/l)	age_y	200	NA	-0.0498634	0.0234825	28.44814	-2.1234293	0.04254	0.05200	.
Delta antibiotic tt days	Baseline (Intercept)	35	236	12.6828934	12.6916885	44.99841	0.9993070	0.32299	0.32299	ns
Delta antibiotic tt days	3	29	236	-11.3029382	7.3110712	197.07435	-1.5460030	0.12371	0.13608	ns
Delta antibiotic tt days	6	30	236	-16.4788730	7.2557827	197.94227	-2.2711365	0.02421	0.03330	*
Delta antibiotic tt days	9	31	236	-29.4923614	7.1817520	197.64750	-4.1065692	0.00006	0.00011	***
Delta antibiotic tt days	12	25	236	-43.9172016	7.7131424	199.65497	-5.6938145	0.00000	0.00000	***
Delta antibiotic tt days	15	23	236	-42.1089966	7.9473336	200.16922	-5.2985062	0.00000	0.00000	***
Delta antibiotic tt days	18	22	236	-46.5567583	8.0877718	201.12244	-5.7564382	0.00000	0.00000	***
Delta antibiotic tt days	21	22	236	-43.8305236	8.0833203	201.42531	-5.4223416	0.00000	0.00000	***
Delta antibiotic tt days	24	19	236	-47.3194488	8.5588867	202.58688	-5.5286921	0.00000	0.00000	***
Delta antibiotic tt days	sex	236	NA	23.4806043	10.8309936	34.14058	2.1679086	0.03723	0.04550	*
Delta antibiotic tt days	age_y	236	NA	-1.0017944	0.3929336	34.92744	-2.5495261	0.01533	0.02410	*

#write_csv(lm_stats_1, "/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/table_outputs/lmm_delta_clinical_par.csv")

2 combine plots for Figure 1

phylum_bars <- readRDS("/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/figures/samples_phylum_wControls.rds")

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

# Arrange the ggplot objects in a grid
grid_arranged <- grid.arrange(
  phylum_bars + labs(tag = "A") + theme(legend.margin = margin(0, 0, 0, 0), legend.box.spacing = unit(0, units = "pt"), legend.box.margin = margin(0, 0, 0, 0)),
  Dchloride + labs(tag = "B"),
  Dppfev1 + labs(tag = "C"),
  Digg + labs(tag = "D"),
  Datb+ labs(tag = "E"),
  ncol = 4, nrow = 2,
  heights = c(1.5, 1),
  widths=c(1.15,1.5,1.5,1.5),
  layout_matrix = rbind(c(1,1,1,1),
                        c(2,3,4,5))
)

# Print or save the arranged plot
grid_arranged

## TableGrob (2 x 4) "arrange": 5 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-4) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
## 3 3 (2-2,2-2) arrange gtable[layout]
## 4 4 (2-2,3-3) arrange gtable[layout]
## 5 5 (2-2,4-4) arrange gtable[layout]

#ggsave("/Users/rebecca/Documents/Forschung/IMMProveCF/Manuscript_16S/Submission_CHM/figure1.tiff",grid_arranged, dpi = 600, width = 17, height = 10)

3 Plot evolution of clinical parameters and run linear mixed effect models for statistical testing

#### Sweat chloride ####
lm <- summary(lmerTest::lmer(chlorid_mmol_cl_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(chlorid_mmol_cl_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_chloride <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_chloride $N_obs <- replace_na(lmm_chloride $N_obs, lmm_chloride $Total_N_obs[1])

lmm_chloride  <- lmm_chloride %>% 
  mutate(Parameter = "Sweat chloride (mmol/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Chloride<- full_metadata_distinct%>%
    ggplot(aes(y=chlorid_mmol_cl_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.5))+
    theme_classic()+
    theme_classic()+
    ylab("Sweat chloride (mmol/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 90, step.increase = 0.05)
Chloride

### ppfev1 ###

lm <- summary(lmerTest::lmer(pp_fev_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(pp_fev_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_ppfev <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_ppfev $N_obs <- replace_na(lmm_ppfev $N_obs, lmm_ppfev $Total_N_obs[1])

lmm_ppfev  <- lmm_ppfev  %>% 
  mutate(Parameter = "ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

ppfev1<- full_metadata_distinct%>%
    ggplot(aes(y=pp_fev_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("ppFEV1")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 145, step.increase = 0.05)
ppfev1

### ppfev1 ###

lm <- summary(lmerTest::lmer(pp_fvc_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(pp_fvc_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_ppfvc  <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_ppfvc$N_obs <- replace_na(lmm_ppfvc$N_obs, lmm_ppfvc$Total_N_obs[1])

lmm_ppfvc  <- lmm_ppfvc   %>% 
  mutate(Parameter = "ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

ppfvc<- full_metadata_distinct%>%
    ggplot(aes(y=pp_fvc_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("ppFVC")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 145, step.increase = 0.05)
ppfvc

#### IgG ####
lm <- summary(lmerTest::lmer(ig_g_g_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(ig_g_g_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_igg <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_igg$N_obs <- replace_na(lmm_igg$N_obs, lmm_igg$Total_N_obs[1])

lmm_igg <- lmm_igg%>% 
  mutate(Parameter = "IgG (g/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

igG <- full_metadata_distinct%>%
    ggplot(aes(y=ig_g_g_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("IgG (g/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,100,200))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 25, step.increase = 0.05)
igG

#### Interleukin 6 ####

lm <- summary(lmerTest::lmer(interleukin6_pg_ml~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(interleukin6_pg_ml)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_il6 <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_il6$N_obs <- replace_na(lmm_il6$N_obs, lmm_il6$Total_N_obs[1])

lmm_il6  <- lmm_il6  %>% 
  mutate(Parameter = "Interleukin 6 (pg/ml)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Il6 <- full_metadata_distinct%>%
    ggplot(aes(y=interleukin6_pg_ml, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Interleukin 6 (pg/ml)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
Il6

#### Interleukin 8 ####

lm <- summary(lmerTest::lmer(interleukin8_pg_ml~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(interleukin8_pg_ml)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_il8 <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_il8$N_obs <- replace_na(lmm_il8$N_obs, lmm_il8$Total_N_obs[1])

lmm_il8  <- lmm_il8  %>% 
  mutate(Parameter = "Interleukin 8 (pg/ml)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

Il8 <- full_metadata_distinct%>%
    ggplot(aes(y=interleukin8_pg_ml, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Interleukin 8 (pg/ml)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
Il8

#### CRP ####

lm <- summary(lmerTest::lmer(crp_mg_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(crp_mg_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_crp <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_crp$N_obs <- replace_na(lmm_crp$N_obs, lmm_crp$Total_N_obs[1])

lmm_crp  <- lmm_crp  %>% 
  mutate(Parameter = "CRP (mg/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

CRP <- full_metadata_distinct%>%
    ggplot(aes(y=crp_mg_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("CRP (mg/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
CRP

#### Calprotectin in stool ####

lm <- summary(lmerTest::lmer(calprotectin_amg_g_stuhl~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(calprotectin_amg_g_stuhl)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_calp <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_calp$N_obs <- replace_na(lmm_calp$N_obs, lmm_calp$Total_N_obs[1])

lmm_calp <- lmm_calp %>% 
  mutate(Parameter = "Calprotectin (mg/g)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

calp_stool <- full_metadata_distinct%>%
    ggplot(aes(y=calprotectin_amg_g_stuhl, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Calprotectin in stool (mg/g)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 3, step.increase = 0.05)
calp_stool

#### GOT/ASAT liver enzyme ####
lm <- summary(lmerTest::lmer(got_asat_u_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(got_asat_u_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_got <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_got$N_obs <- replace_na(lmm_got$N_obs, lmm_got$Total_N_obs[1])

lmm_got <- lmm_got %>% 
  mutate(Parameter = "GOT/ASAT (U/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

got<- full_metadata_distinct%>%
    ggplot(aes(y=got_asat_u_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("GOT/ASAT (ul/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2, step.increase = 0.05)
got

#### GPT/ALAT liver enzyme ####

lm <- summary(lmerTest::lmer(gpt_alat_u_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(gpt_alat_u_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_gpt <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_gpt$N_obs <- replace_na(lmm_gpt$N_obs, lmm_gpt$Total_N_obs[1])

lmm_gpt <- lmm_gpt %>% 
  mutate(Parameter = "GPT/ALAT (U/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

gpt<- full_metadata_distinct%>%
    ggplot(aes(y=gpt_alat_u_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("GPT/ALAT (ul/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2, step.increase = 0.025)
gpt

### HbA1c ###

lm <- summary(lmerTest::lmer(hb_a1c_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(hb_a1c_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_hba1c <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_hba1c$N_obs <- replace_na(lmm_hba1c$N_obs, lmm_hba1c$Total_N_obs[1])

lmm_hba1c <- lmm_hba1c %>% 
  mutate(Parameter = "HbA1c (%)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

hba1c<- full_metadata_distinct%>%
    ggplot(aes(y=hb_a1c_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("HbA1c (%)")+
    xlab("Months from ETI treatment start")+
    scale_y_log10()+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 1.1, step.increase = 0.025)
hba1c

### Leucocytes ###

lm <- summary(lmerTest::lmer(leukozyten_nl~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))

## New names:
## • `` -> `...6`

lm_stats

# print lm stats into data frame and add N of observations

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(leukozyten_nl)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24"))%>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_leukos <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_leukos$N_obs <- replace_na(lmm_leukos$N_obs, lmm_leukos$Total_N_obs[1])

lmm_leukos <- lmm_leukos %>% 
  mutate(Parameter = "Leukocytes (/nl)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))

## New names:
## • `` -> `...9`

leucos<- full_metadata_distinct%>%
    ggplot(aes(y=leukozyten_nl, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Leucocytes (/nl)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 1.4, step.increase = 0.05)
leucos

ggarrange(ppfev1, ppfvc, Chloride, igG, Il6, calp_stool, gpt, hba1c, nrow = 2, ncol = 4)

#ggsave("/Users/rebecca/Documents/Forschung/IMMProveCF/Manuscript_16S/Submission_CHM/supp_figure1.tiff", width = 14, height = 7)

stats_table <- rbind(lmm_chloride,lmm_Dppfev,lmm_ppfvc,lmm_igg,lmm_leukos, lmm_il6, lmm_il8, lmm_calp, lmm_gpt, lmm_got, lmm_hba1c)
#write_csv(stats_table, "/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/table_outputs/lmm_clinical_parameters.csv")
kable(stats_table)

Parameter	Months_after_ETI_start	N_obs	Total_N_obs	Estimate	Std..Error	df	t.value	p	fdr	fdr_star
Sweat chloride (mmol/l)	Baseline (Intercept)	34	119	69.0854861	4.3279602	42.47775	15.9625974	0.00000	0.00000	***
Sweat chloride (mmol/l)	3	25	119	-35.6683953	2.7298654	82.09485	-13.0659903	0.00000	0.00000	***
Sweat chloride (mmol/l)	6	6	119	-35.0235154	4.9688991	91.98396	-7.0485463	0.00000	0.00000	***
Sweat chloride (mmol/l)	9	3	119	-37.3910610	6.8785114	91.54992	-5.4359234	0.00000	0.00000	***
Sweat chloride (mmol/l)	12	19	119	-40.6820780	3.0234061	84.78035	-13.4557106	0.00000	0.00000	***
Sweat chloride (mmol/l)	15	5	119	-39.2537245	5.3309544	90.11204	-7.3633578	0.00000	0.00000	***
Sweat chloride (mmol/l)	18	14	119	-44.5773069	3.4142706	85.44003	-13.0561726	0.00000	0.00000	***
Sweat chloride (mmol/l)	21	2	119	-32.9562276	8.1513534	88.86314	-4.0430375	0.00011	0.00014	***
Sweat chloride (mmol/l)	24	11	119	-42.0315021	3.7597756	86.30286	-11.1792582	0.00000	0.00000	***
Sweat chloride (mmol/l)	sex	119	NA	10.0339457	3.6951540	32.75024	2.7154337	0.01048	0.01048	*
Sweat chloride (mmol/l)	age_y	119	NA	0.4255756	0.1348529	33.68501	3.1558499	0.00336	0.00370	**
Delta ppFEV1	Baseline (Intercept)	35	235	1.0158422	3.4045717	38.64296	0.2983759	0.76701	0.84372	ns
Delta ppFEV1	3	30	235	12.4628627	1.5482261	193.34212	8.0497692	0.00000	0.00000	***
Delta ppFEV1	6	31	235	11.3119118	1.5384045	194.31464	7.3530151	0.00000	0.00000	***
Delta ppFEV1	9	31	235	13.4351520	1.5389186	194.47907	8.7302550	0.00000	0.00000	***
Delta ppFEV1	12	26	235	9.3934309	1.6321052	196.03302	5.7554077	0.00000	0.00000	***
Delta ppFEV1	15	22	235	10.2279562	1.7298648	196.79571	5.9125754	0.00000	0.00000	***
Delta ppFEV1	18	20	235	7.4812868	1.7930864	198.43422	4.1722956	0.00005	0.00008	***
Delta ppFEV1	21	22	235	6.8596951	1.7354919	198.80899	3.9525942	0.00011	0.00017	***
Delta ppFEV1	24	18	235	6.3404820	1.8748622	200.42886	3.3818389	0.00087	0.00119	**
Delta ppFEV1	sex	235	NA	-1.4902106	2.9866734	32.30319	-0.4989533	0.62119	0.75924	ns
Delta ppFEV1	age_y	235	NA	-0.0100933	0.1079618	33.07653	-0.0934894	0.92608	0.92608	ns
ppFEV1	Baseline (Intercept)	35	234	101.1293680	5.8076667	34.29655	17.4130806	0.00000	0.00000	***
ppFEV1	3	30	234	10.6705013	1.5983305	190.94312	6.6760293	0.00000	0.00000	***
ppFEV1	6	31	234	9.5525085	1.5907892	192.46793	6.0048864	0.00000	0.00000	***
ppFEV1	9	31	234	11.6134594	1.5922728	193.23643	7.2936368	0.00000	0.00000	***
ppFEV1	12	26	234	9.4728823	1.6929252	195.44954	5.5955706	0.00000	0.00000	***
ppFEV1	15	21	234	8.2862901	1.8252527	197.22898	4.5398043	0.00001	0.00002	***
ppFEV1	18	20	234	6.7752830	1.8704931	201.28564	3.6221908	0.00037	0.00051	***
ppFEV1	21	22	234	6.7483801	1.8122199	202.30180	3.7238197	0.00025	0.00040	***
ppFEV1	24	18	234	4.9188976	1.9654334	205.93999	2.5027038	0.01310	0.01601	*
ppFEV1	sex	234	NA	5.4206413	5.2596917	31.46821	1.0306006	0.31058	0.31058	ns
ppFEV1	age_y	234	NA	-0.3240456	0.1884772	33.01347	-1.7192824	0.09493	0.10442	ns
IgG (g/l)	Baseline (Intercept)	30	232	10.2534962	1.2180677	34.63855	8.4178378	0.00000	0.00000	***
IgG (g/l)	3	29	232	-1.2524893	0.2837596	189.26840	-4.4139106	0.00002	0.00002	***
IgG (g/l)	6	31	232	-1.7204948	0.2818769	191.18140	-6.1037089	0.00000	0.00000	***
IgG (g/l)	9	31	232	-1.8475817	0.2825581	192.28914	-6.5387675	0.00000	0.00000	***
IgG (g/l)	12	26	232	-1.9040301	0.3005890	195.01382	-6.3343315	0.00000	0.00000	***
IgG (g/l)	15	22	232	-1.9622117	0.3183205	197.82439	-6.1642640	0.00000	0.00000	***
IgG (g/l)	18	22	232	-1.7266318	0.3208876	202.32241	-5.3807989	0.00000	0.00000	***
IgG (g/l)	21	22	232	-1.8791860	0.3214042	203.76324	-5.8467995	0.00000	0.00000	***
IgG (g/l)	24	19	232	-1.5403600	0.3402355	208.91737	-4.5273346	0.00001	0.00001	***
IgG (g/l)	sex	232	NA	-1.3508350	1.1103225	31.91779	-1.2166150	0.23267	0.23267	ns
IgG (g/l)	age_y	232	NA	0.1321313	0.0395804	34.22812	3.3383045	0.00204	0.00225	**
Leukocytes (/nl)	Baseline (Intercept)	32	234	9.3385414	1.2322962	73.25221	7.5781629	0.00000	0.00000	***
Leukocytes (/nl)	3	30	234	-1.8968711	1.1078251	199.40245	-1.7122478	0.08841	0.20730	ns
Leukocytes (/nl)	6	30	234	-2.0856793	1.1105841	201.01311	-1.8780021	0.06183	0.20730	ns
Leukocytes (/nl)	9	30	234	-1.7918893	1.1099385	200.76492	-1.6144042	0.10801	0.20730	ns
Leukocytes (/nl)	12	26	234	-1.7316513	1.1595152	203.45800	-1.4934270	0.13687	0.20730	ns
Leukocytes (/nl)	15	23	234	-1.6622048	1.2059396	203.11072	-1.3783484	0.16961	0.20730	ns
Leukocytes (/nl)	18	22	234	0.6454324	1.2252530	204.03064	0.5267748	0.59892	0.65881	ns
Leukocytes (/nl)	21	22	234	-1.7581928	1.2239677	204.27608	-1.4364699	0.15240	0.20730	ns
Leukocytes (/nl)	24	19	234	-1.8208489	1.2918426	205.88786	-1.4094975	0.16020	0.20730	ns
Leukocytes (/nl)	sex	234	NA	-1.7000438	0.9172221	34.43945	-1.8534701	0.07240	0.20730	ns
Leukocytes (/nl)	age_y	234	NA	-0.0004959	0.0336626	36.37515	-0.0147317	0.98833	0.98833	ns
Interleukin 6 (pg/ml)	Baseline (Intercept)	23	142	12.0913248	3.3496209	66.71004	3.6097591	0.00059	0.00527	**
Interleukin 6 (pg/ml)	3	17	142	-10.4236933	3.1643478	119.79613	-3.2941048	0.00130	0.00527	**
Interleukin 6 (pg/ml)	6	14	142	-8.2930616	3.3910124	123.08216	-2.4456005	0.01588	0.03283	*
Interleukin 6 (pg/ml)	9	19	142	-10.0901876	3.0916936	118.85769	-3.2636442	0.00144	0.00527	**
Interleukin 6 (pg/ml)	12	16	142	-7.7979661	3.2482240	119.60932	-2.4006861	0.01791	0.03283	*
Interleukin 6 (pg/ml)	15	7	142	-8.3825303	4.3196965	120.45984	-1.9405368	0.05465	0.06876	.
Interleukin 6 (pg/ml)	18	12	142	-6.6609438	3.5554994	120.25715	-1.8734200	0.06344	0.06978	.
Interleukin 6 (pg/ml)	21	19	142	-5.9474023	3.0852043	119.86489	-1.9277175	0.05626	0.06876	.
Interleukin 6 (pg/ml)	24	15	142	-9.9125157	3.3380392	118.40559	-2.9695624	0.00361	0.00993	**
Interleukin 6 (pg/ml)	sex	142	NA	-5.5781142	2.4148587	31.84047	-2.3099133	0.02753	0.04326	*
Interleukin 6 (pg/ml)	age_y	142	NA	0.1478922	0.0894327	38.45241	1.6536708	0.10634	0.10634	ns
Interleukin 8 (pg/ml)	Baseline (Intercept)	28	212	33.4005300	12.6078088	101.07647	2.6491939	0.00936	0.10301	ns
Interleukin 8 (pg/ml)	3	24	212	-13.6132223	13.6578198	183.33274	-0.9967347	0.32021	0.62451	ns
Interleukin 8 (pg/ml)	6	29	212	-13.4123485	13.0144032	183.94510	-1.0305773	0.30409	0.62451	ns
Interleukin 8 (pg/ml)	9	28	212	-10.1635654	13.1531518	185.83268	-0.7727095	0.44068	0.62451	ns
Interleukin 8 (pg/ml)	12	23	212	-0.8228657	13.9117605	188.35988	-0.0591489	0.95290	0.95290	ns
Interleukin 8 (pg/ml)	15	21	212	12.5951739	14.3039541	188.67320	0.8805379	0.37969	0.62451	ns
Interleukin 8 (pg/ml)	18	21	212	-1.2850447	14.3494815	189.41538	-0.0895534	0.92874	0.95290	ns
Interleukin 8 (pg/ml)	21	22	212	10.5982778	14.1311375	189.24140	0.7499947	0.45419	0.62451	ns
Interleukin 8 (pg/ml)	24	16	212	-12.9644453	15.6343884	191.24249	-0.8292263	0.40801	0.62451	ns
Interleukin 8 (pg/ml)	sex	212	NA	7.4896144	8.3854131	33.85495	0.8931718	0.37807	0.62451	ns
Interleukin 8 (pg/ml)	age_y	212	NA	-0.1625423	0.3145589	37.69665	-0.5167308	0.60836	0.74355	ns
Calprotectin (mg/g)	Baseline (Intercept)	18	120	178.7703785	25.9977966	50.57516	6.8763665	0.00000	0.00000	***
Calprotectin (mg/g)	3	9	120	-117.3008502	34.9875511	104.56807	-3.3526453	0.00112	0.00153	**
Calprotectin (mg/g)	6	14	120	-106.2924317	30.2901786	97.20218	-3.5091385	0.00068	0.00107	**
Calprotectin (mg/g)	9	15	120	-88.1664265	30.0300230	92.32422	-2.9359427	0.00420	0.00513	**
Calprotectin (mg/g)	12	11	120	-126.6824229	32.8491393	97.59591	-3.8564914	0.00021	0.00038	***
Calprotectin (mg/g)	15	11	120	-139.7518941	33.0813239	97.17936	-4.2244952	0.00005	0.00012	***
Calprotectin (mg/g)	18	12	120	-143.6487962	32.1672533	105.83157	-4.4656842	0.00002	0.00006	***
Calprotectin (mg/g)	21	16	120	-159.0050397	29.5171626	99.20808	-5.3868674	0.00000	0.00000	***
Calprotectin (mg/g)	24	14	120	-164.8620305	30.9167393	100.48976	-5.3324521	0.00000	0.00000	***
Calprotectin (mg/g)	sex	120	NA	-12.0107707	18.4637036	19.30519	-0.6505071	0.52303	0.52303	ns
Calprotectin (mg/g)	age_y	120	NA	0.9174243	0.6591936	17.76059	1.3917374	0.18119	0.19931	ns
GPT/ALAT (U/l)	Baseline (Intercept)	33	236	12.8189870	5.6236280	44.68613	2.2794870	0.02747	0.04316	*
GPT/ALAT (U/l)	3	30	236	11.6287612	3.4478332	195.46765	3.3727737	0.00090	0.00445	**
GPT/ALAT (U/l)	6	31	236	10.7105478	3.4252042	196.50848	3.1269809	0.00203	0.00559	**
GPT/ALAT (U/l)	9	31	236	7.9681772	3.4261340	196.56343	2.3257051	0.02105	0.03860	*
GPT/ALAT (U/l)	12	26	236	11.9048006	3.6264203	198.53358	3.2827967	0.00121	0.00445	**
GPT/ALAT (U/l)	15	23	236	9.4939656	3.7839403	199.08265	2.5090157	0.01290	0.02839	*
GPT/ALAT (U/l)	18	21	236	6.4525604	3.9080376	200.04581	1.6510999	0.10029	0.13789	ns
GPT/ALAT (U/l)	21	22	236	3.2229866	3.8473249	200.33394	0.8377215	0.40319	0.40319	ns
GPT/ALAT (U/l)	24	19	236	6.2267434	4.0710685	201.54846	1.5295108	0.12771	0.15608	ns
GPT/ALAT (U/l)	sex	236	NA	18.2526070	4.7434809	32.26636	3.8479352	0.00053	0.00445	**
GPT/ALAT (U/l)	age_y	236	NA	0.1657523	0.1723532	33.22639	0.9617012	0.34315	0.37746	ns
GOT/ASAT (U/l)	Baseline (Intercept)	33	236	27.0033788	2.9279484	42.07908	9.2226281	0.00000	0.00000	***
GOT/ASAT (U/l)	3	30	236	3.1357365	1.6487210	194.82588	1.9019207	0.05866	0.21508	ns
GOT/ASAT (U/l)	6	31	236	2.7475023	1.6385718	195.89342	1.6767665	0.09518	0.26175	ns
GOT/ASAT (U/l)	9	31	236	1.6284189	1.6390431	195.97009	0.9935180	0.32168	0.58975	ns
GOT/ASAT (U/l)	12	26	236	1.7633046	1.7357418	197.80044	1.0158795	0.31093	0.58975	ns
GOT/ASAT (U/l)	15	22	236	1.4495953	1.8381258	198.49184	0.7886268	0.43127	0.59300	ns
GOT/ASAT (U/l)	18	22	236	0.5760887	1.8435197	199.47163	0.3124939	0.75499	0.86463	ns
GOT/ASAT (U/l)	21	22	236	0.5008440	1.8424279	199.80357	0.2718391	0.78603	0.86463	ns
GOT/ASAT (U/l)	24	19	236	-0.0988094	1.9502776	201.11846	-0.0506643	0.95964	0.95964	ns
GOT/ASAT (U/l)	sex	236	NA	9.7185740	2.5047824	32.00660	3.8800072	0.00049	0.00269	**
GOT/ASAT (U/l)	age_y	236	NA	-0.0800732	0.0909004	32.90946	-0.8808894	0.38477	0.59300	ns
HbA1c (%)	Baseline (Intercept)	16	123	5.9435974	0.5022168	39.10772	11.8347242	0.00000	0.00000	***
HbA1c (%)	3	19	123	-0.6201108	0.1799150	83.35358	-3.4466875	0.00089	0.00490	**
HbA1c (%)	6	12	123	-0.6172339	0.2286378	87.89141	-2.6996145	0.00833	0.03053	*
HbA1c (%)	9	12	123	-0.5039793	0.2060707	83.87834	-2.4456628	0.01655	0.03641	*
HbA1c (%)	12	11	123	-0.5162835	0.2092032	83.60846	-2.4678569	0.01563	0.03641	*
HbA1c (%)	15	9	123	-0.4156051	0.2268026	83.73601	-1.8324528	0.07044	0.11069	ns
HbA1c (%)	18	10	123	-0.4834953	0.2131365	84.69849	-2.2684777	0.02584	0.04738	*
HbA1c (%)	21	19	123	-0.2903652	0.1779156	85.67842	-1.6320387	0.10634	0.14622	ns
HbA1c (%)	24	15	123	-0.0533455	0.1862931	86.86123	-0.2863524	0.77529	0.77529	ns
HbA1c (%)	sex	123	NA	-0.3215008	0.4365231	32.20279	-0.7365035	0.46676	0.51343	ns
HbA1c (%)	age_y	123	NA	0.0231887	0.0158975	33.78674	1.4586412	0.15390	0.18810	ns